program d04(input,output);

uses common,d04_header;
uses relax,tools,etools,inp,outp,ppot,decor,tile_mc,diffract;

{Project started Pittsburgh-Ithaca, summer 2000}
{This revision Ithaca 2002}

{
### INPUT FILES:
d04.inp
uc.inp
pp.inp
objects.def
[atoms.def]
[orbit_labels.def]
### OUTPUT FILES:
copa.out
copa.min
copas.out
tiling.out
tiling.min
tilings.out
objects.ini
objects.out
hc.rep
u.rep
...u
}
Const
   nnbond = 10; {number of bond types}
   nnln	  = 800; {number of site neighbors within the cutoff dist}
   nnpair = nsites*70; {pairs for swaps}
   nlh	  = nnln div 2; {actually used for nne, since i-j are symmetric}
Type
   rdtt	   = array[1..ndtt] of myreal; {ndtt comes from ppot: number of pairs.}
   rbonds  = array[1..nnbond] of v2;
   chbonds = array[1..nnbond] of i2;
   tltl	   = array[1..t_ele,1..t_ele] of myreal;
Var
   pl	       : array[1..nsites] of i2;
   nnn,s_sites : array[1..nsites] of integer;
   s_frozen    : array[1..nsites] of boolean;
   nnl,nl2     : array[1..nsites,1..nnln] of integer;
   nne	       : array[1..nsites*nlh] of tltl; {assuming at most ternary}
   lpair       : array[1..nnpair,1..2] of integer;
   npp	       : array[1..nnpair] of integer;
   {rpair      : array[1..nnpair] of myreal;}
   bpair       : array[1..nnpair,1..nnbond] of boolean;
Var    {misc declarations}   
   nfree,nij,imcs	     : integer; {number of unfrozen sites}
   loops,hiter,viter,relmcs0 : integer;
   hiterw		     : integer;
   relmcs,lmcs,smcs,smcs2    : integer;
   loops_ini,gmcs,tmcs,omcs  : integer;
   twopa		     : array[1..10] of i2;
   dum			     : myreal;
   chr			     : char;
   eofile		     : boolean;
   b			     : i5;
   u,v,w		     : v2;
   txt,txi,txyzs,txef	     : text;
   txpardf,txrepdf	     : text;
   i,j,k,l,m,n,o,p	     : integer; 
   code,test_glob,test_perm  : integer;
   nswapm		     : integer;
   no1,at1,no2,at2	     : array[1..200] of integer;
   atr,ato1,ato2	     : array[1..200] of integer;
   e0,e1		     : myreal;
   cmp,por		     : MYREAL;		       
   x2			     : i2;				       
   lgchp,eval_nom,rchemi     : boolean;
   rchemx,cheminp	     : boolean;
   flip_nu		     : array[1..20] of i2;		       
   nflipnu,level,flipmode    : integer;			       
   xbeta		     : myreal;
   rcuts,sigmc		     : myreal;
   nbond,nbpa,tlevel	     : integer;
   bondch		     : chbonds;
   bondr		     : rbonds;
   couple,vpin		     : integer;
   control,control2	     : integer;
   xmele		     : array[1..t_ele] of myreal;
   kboltz,ecur		     : myreal;
   nmax			     : integer;
   rcutpair,vol,rho,rho_lim  : myreal;
   {NELE is defined in ppot (=z_ele), z_ele in common.p (=110)}
   u_self		     : array[0..nele] of myreal;
   chpt			     : array[0..nele] of myreal;
   npair		     : integer;
   nchs,natnom,natlim	     : integer;
   chlist		     : array[1..nsites] of integer; {ch-buffer}
   write_stout,table_reinit  : boolean;
   osta			     : array[1..t_obj] of integer;
   fene_tot		     : myreal;
   xgamma,dramax,ueff	     : myreal;
   nmrel,nmrelf		     : integer;
   w3			     : v3;
   ticoname,inpname	     : mystring;
   tinp,nschemes	     : integer;
   sch			     : array[1..100] of i2;
   {DFIT}
   itmp0,dmcsteps        : integer;
   se,sel:rsites;
   fobj:integer;
   io,jo:integer; {misc}
   elas:myreal;
   savesigma:boolean;

procedure site_occ;
Var
   i,l1,l2 : integer;
begin
   ntocc:=ntocc+1;
   for i:=1 to sites do if (ich[i]>0) then begin
      l1:=pl[i,1]; l2:=pl[i,2]; p:=pele[ich[i]];
      n_oc[l1,l2,p]:=n_oc[l1,l2,p]+1;
   end;
end; { site_oc }
procedure list2atoms;
Var
   i,l1,l2 : integer;
begin
   for i:=1 to sites do begin
      l1:=pl[i,1]; l2:=pl[i,2];
      n_a[l1,l2,1]:=ich[i];
      para2xyz(n_para[l1,l2],ich[i],1,1,n_r[l1,l2]);
   end;
end; { list2atoms }
procedure ch2ch(i,j: integer; Var ch1,ch2:integer);
begin
   if (i<j) then begin ch1:=ich[i]; ch2:=ich[j] end else
   begin ch2:=ich[i]; ch1:=ich[j] end;
end;
procedure site_ene_list(po,nn: integer; po10:i10;Var u,up:myreal);
{avoid double-counting}
var
   i,j,o,pe : integer;
   ch1,ch2  : integer;
   ep	    : myreal;
begin
   u:=0.0;
   if (ich[po]>0) then for i:=1 to nnn[po] do begin
      if (ich[nnl[po,i]]>0) then begin
	 pe:=nl2[po,i];
	 ch2ch(po,nnl[po,i],ch1,ch2);
	 ep:=nne[pe,pele[ch1],pele[ch2]];
	 u:=u+ep;
	 for j:=1 to nn do if (nnl[po,i]=po10[j]) then begin
	    up:=up+ep;
	 end;
      end;
   end;
end; { site_ene }
procedure site_ene(po : integer; Var u:myreal);
var
   i,ch1,ch2  : integer;
   p,pn,p2,pe : integer;
   b1,b2      : boolean;
begin
   u:=0.0; if (ich[po]>0) then for i:=1 to nnn[po] do begin
      if (ich[nnl[po,i]]>0) then begin
	 p2:=nnl[po,i];
	 ch2ch(po,p2,ch1,ch2);
	 pe:=nl2[po,i];
	 u:=u+nne[pe,pele[ch1],pele[ch2]];
      end;
   end;
end; { site_ene }
function rin(p :integer; rr:myreal):boolean;
Var
   rin0	 : boolean;
begin
   rin0:=false;
   if ((rr>bondr[p,1]) and (rr<bondr[p,2])) then rin0:=true;
   rin:=rin0;
end; { rin }
procedure atom_list;
Var
   i,j,o : integer;
begin
   sites:=0;
   for i:=1 to tnodes do begin
      for j:=1 to nan[i] do begin
	 sites:=sites+1; n_a[i,j,3]:=sites;
	 if (sites>nsites) then begin
	    writeln('Insufficient array dimension for site-list.');
	    writeln('Current Const nsites = ',nsites:1);
	    halt;
	 end;
	 if o_frozen[n_a[i,j,2]] then
	    s_frozen[sites]:=true else s_frozen[sites]:=false;
	 pl[sites,1]:=i; pl[sites,2]:=j;
	 och[sites]:=n_a[i,j,2];
	 ich[sites]:=n_a[i,j,1];
	 lay[sites]:=((i-1) div tnodes0);
      end;
   end;
end; { atom_list }
procedure enn_table;
{Chemical displacements break the R(i,j)==R(j,i) symmetry!!!
(Even when atoms belong to the same orbit, the tile-objects dont
have necessarily the same orientation!)
We store U(R(A,B)) in the upper triangle, and U(R(B,A)) in the lower.
}
Const
   iucp	 = 4;
   iucpp = (iucp*2+1)*(iucp*2+1)*(iucp*2+1); {for flat lists gen. by 3 iext loops}
Var
   i,j,l,o,p1,p2,p3 : integer;
   orb,o1,o2	    : integer;
   l1,l2,npx,k1,k2  : integer;
   ar,u,ud,udd,xx   : myreal;
   eps,ep1,rcpl	    : myreal;
   v,r1,r2,w,va,vex : v3;
   new_pair	    : array[1..iucpp] of boolean;
   px		    : array[1..iucpp] of v3;
   ide,first	    : boolean;
   dra		    : array[1..nsites,1..t_ele] of v3;
   ara		    : array[1..t_ele,1..t_ele] of myreal;
   sr		    : array[1..nsites] of v3;
   p		    : array[1..t_orb,1..z_ele] of integer;
   ofa		    : array[1..t_ele,1..t_ele] of boolean;
begin
   ATOM_LIST;
   eps:=1e-6;
   rcpl:=0.0; {in Angstroms}
   for i:=1 to t_orb do
      for j:=1 to z_ele do p[i,j]:=0;
   for i:=1 to sites do begin {precalculate chemical displacements}
      orb:=och[i];
      l1:=pl[i,1]; l2:=pl[i,2];
      para2xyz(n_para[l1,l2],ich[i],1,0,v);
      for o:=1 to 3 do sr[i,o]:=bound(v[o]+rt[l1,o]);
      for j:=1 to noche[orb] do begin
	 para2xyz(n_para[l1,l2],oche[orb,j],0,1,v);
	 ortho3(v,dra[i,j]); {!! in [A] !!}
	 xx:=absv3(dra[i,j]);
	 if (xx>eps) then p[orb,oche[orb,j]]:=j;
	 if (xx>rcpl/2.0) then rcpl:=2.0*xx; {extend cutoff by this}
      end;
      if debug[4] then begin
	 para2xyz(n_para[l1,l2],ich[i],1,1,v);
	 if (dist3(n_r[l1,l2],v)>1e-7) then begin
	    write('Site ',l1:1,' ',l2:1,': ');
	    writeln(' wrong cur.pos. n_r ');
	    halt;
	 end;
      end;
   end;
   npx:=0;
   for o:=1 to nele do u_self[o]:=0.0; {self-energy for LG}
   for p1:=-iext[1] to iext[1] do
      for p2:=-iext[2] to iext[2] do
	 for p3:=-iext[3] to iext[3] do begin
	    ide:=((p1=0) and (p2=0) and (p3=0));
	    npx:=npx+1;
	    if(npx>iucpp)then begin
	       writeln('ENN-table in d04_.p: need to increase const iucpp.');
	       halt;
	    end;
	    new_pair[npx]:=ide;
	    px[npx,1]:=p1; px[npx,2]:=p2; px[npx,3]:=p3;
	    ortho3(px[npx],w);
	    px[npx]:=w; {!! in [A] !!}
	    ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
	    if ((ar<rcut) and not(ide)) then begin
	       for o:=1 to nnele do begin
		  upp(ele[o],ele[o],ar,u,ud,udd);
		  u_self[ele[o]]:=u_self[ele[o]]+u/2.0;
	       end;
	    end;
	 end;
   npair:=0; nij:=0;
   for o:=1 to nnbond do npp[o]:=0;
   nmax:=0;
   for i:=1 to nsites*nlh do
      for p1:=1 to t_ele do
	 for p2:=1 to t_ele do nne[i,p1,p2]:=0.0;
   for i:=1 to nnpair do
      for j:=1 to nnbond do bpair[i,j]:=false;
   for i:=1 to sites do nnn[i]:=0;
   for i:=1 to sites-1 do begin
      for j:=i+1 to sites do begin
	 for o:=1 to 3 do v[o]:=bound(sr[j,o]-sr[i,o]);
	 ortho3(v,va);
	 first:=true;
	 for p1:=1 to npx do begin
	    for o:=1 to 3 do vex[o]:=va[o]+px[p1,o];
	    ar:=sqrt(vex[1]*vex[1]+vex[2]*vex[2]+vex[3]*vex[3]);
	    if (ar<(rcut+rcpl)) then begin {rcpl: for chem-displ.}
	       if (first) then begin {avoid repeated entries for the same pair}
		  nnn[i]:=nnn[i]+1;
		  nnn[j]:=nnn[j]+1;
		  if ((nnn[i]>nnln) or (nnn[j]>nnln)) then begin
		     writeln('Enn table: insufficient array dimension for neighbor tables. Current nnln=',nnln:1,'.');
		     halt;
		  end;
		  nnl[i,nnn[i]]:=j;
		  nnl[j,nnn[j]]:=i;
		  nij:=nij+1;
		  if (nij>nsites*nlh) then begin
		     writeln('Enn table: insufficient array dim for tabulated pair energies. Increase # of sites NSITES or bonds/site NNLN');
		     halt;
		  end;
		  nl2[i,nnn[i]]:=nij;
		  nl2[j,nnn[j]]:=nij;
		  if (nnn[i]>nmax) then nmax:=nnn[i];
		  if (nnn[j]>nmax) then nmax:=nnn[j];
		  first:=false;
	       end;
	       if new_pair[p1] then begin {SHORTEST BOND FOR A PAIR ONLY!}
		  if not(s_frozen[i] or s_frozen[j]) then begin {tab of pairs}
		     for o:=1 to nbond do if rin(o,ar) then begin
			npair:=npair+1;
			if (npair>nnpair) then begin
			   writeln('npair>nnpair: increase Const nnpair in d04_.p.');
			   halt;
			end;
			lpair[npair,1]:=i;
			lpair[npair,2]:=j;
			bpair[npair,o]:=true;
			npp[o]:=npp[o]+1;
		     end;
		  end;
	       end;
	       {NOTE: with chem-displacement, U(i,j) is no longer symmetric!!}
	       {Note the IF cond. in site_ene for accessing nne[nij,i,j]}
	       for l1:=1 to t_ele do
		  for l2:=1 to t_ele do ara[l1,l2]:=ar;
	       for l1:=1 to nnele do begin
		  for l2:=1 to nnele do begin
		     k1:=p[och[i],ele[l1]];
		     k2:=p[och[j],ele[l2]];
		     if ((k1>0) or (k2>0)) then begin
			for o:=1 to 3 do w[o]:=vex[o];
			if (k1>0) then
			   for o:=1 to 3 do w[o]:=w[o]-dra[i,k1,o];
			if (k2>0) then
			   for o:=1 to 3 do w[o]:=w[o]+dra[j,k2,o];
			ara[l1,l2]:=absv3(w);
		     end;
		     u:=0.0;
		     if(ara[l1,l2]<rcut)then upp(ele[l1],ele[l2],ara[l1,l2],u,ud,udd);
		     nne[nij,l1,l2]:=nne[nij,l1,l2]+u;
		  end;
	       end;
	    end;
	 end;
      end;
   end;
   if debug[4] then begin
      ep1:=0.0; k:=0;
      for i:=1 to sites do begin
	 if (ich[i]>0) then k:=k+1;
	 site_ene(i,xx); ep1:=ep1+u_self[ich[i]]+xx/2.0;
      end;
      ep1:=ep1/k;
      writeln('Enn_table: E=',ep1:1:6,' [eV/atom]');
      writeln('Enn table: Max. ',nmax:1,' neighbs, ave ',nij/k*2:1:3,', ',nij:1,' pairs');
      writeln('Enn_table: rcpl=',rcpl:1:2,'[A]');
   end;
end;
procedure SWAP_two_pairs(kkl,mcs,lab,pb1,pb2 : integer;Var epot,rate:myreal);
{nn-tables, precalc energies...}
Var
   i,j,l,o,nppo	     : integer;
   o1,o2,o3,o4	     : integer;
   l1,l2,k1,k2	     : integer;
   p1,p2,po1,po2     : integer;
   kk		     : integer;
   ar,x,y	     : myreal;
   u0,u1,du0,du	     : myreal;
   v,r1,r2,w	     : v3;
   flptot,nflp	     : integer;
   swappable,identic : boolean;
   check	     : boolean;
   mfact,xp	     : myreal;
   pop1,pop2	     : array[1..nnpair] of integer; {pointers to useful pairs}
   npop1,npop2	     : integer;
   found	     : boolean;
   ppo		     : i10;
   
function chbin(p,ich1,ich2 : integer):boolean;
begin
   chbin:=false;
   if ((ich1=bondch[p,1]) and (ich2=bondch[p,2])) then chbin:=true;
   if ((ich1=bondch[p,2]) and (ich2=bondch[p,1])) then chbin:=true;
   {-1 is wildcard:}
   if ((ich1=bondch[p,1]) and (bondch[p,2]=-1)) then chbin:=true;
   if ((ich1=bondch[p,2]) and (bondch[p,1]=-1)) then chbin:=true;
   if ((ich2=bondch[p,1]) and (bondch[p,2]=-1)) then chbin:=true;
   if ((ich2=bondch[p,2]) and (bondch[p,1]=-1)) then chbin:=true;
   if ((bondch[p,1]=-1) and (bondch[p,2]=-1)) then chbin:=true;
end; { chbin }
begin
   occupied_sites(natoc);
   if debug[5] then writeln('SWAP_two_pairs: E_ini = ',epot/natoc:1:6);
   {list of bonds of type=pb}
   npop1:=0; for i:=1 to npair do if ((bpair[i,pb1]) and (chbin(pb1,ich[lpair[i,1]],ich[lpair[i,2]]))) then begin
      npop1:=npop1+1; pop1[npop1]:=i;
   end;
   npop2:=0; for i:=1 to npair do if ((bpair[i,pb2]) and (chbin(pb2,ich[lpair[i,1]],ich[lpair[i,2]]))) then begin
      npop2:=npop2+1; pop2[npop2]:=i;
   end;
   {for i:=1 to npop1 do begin
   writeln(i:1,' ',ich[lpair[pop1[i],1]]:1,' ',ich[lpair[pop1[i],2]]:1);
   end;
   for i:=1 to npop2 do begin
   writeln(i:1,' ',ich[lpair[pop2[i],1]]:1,' ',ich[lpair[pop2[i],2]]:1);
   end;}
   {in the destroy-recreate-list reshuffling method (as opposed to
    swap-site-list method) we have to demand that only the atoms belonging
    to the same ORBIT may swap; otherwise destruction-reconstruction of a motif
    bound to an object does not conserve chemistry:}
   check:=false; {will check equality of orbit# and layer# for a swap}
   {lets go for it}
   nflp:=0; flptot:=0;
   du:=0.0;
   if debug[5] then writeln('mcs npop1 npop2 nmax : ',mcs:1,' ',npop1:1,' ',npop2:1,' ',nmax:1);
   for kk:=1 to mcs*(npop1+npop2) do begin
      eavnn[kkl,0]:=eavnn[kkl,0]+1.0;
      eavnn[kkl,lab]:=eavnn[kkl,lab]+1.0;
      l:=0; repeat l:=l+1;
	 x:=randu(seed);
	 p1:=trunc(x*npop1)+1;
	 if (p1>npop1) then p1:=npop1;
	 po1:=pop1[p1];
	 l1:=lpair[po1,1];
	 l2:=lpair[po1,2];
	 x:=randu(seed);
	 p2:=trunc(x*npop2)+1;
	 if p2>npop2 then p2:=npop2;
	 po2:=pop2[p2];
	 k1:=lpair[po2,1];
	 k2:=lpair[po2,2];
	 found:=((k1<>l1) and (k1<>l2) and (k2<>l1) and (k2<>l2));
	 {we need to test again chemistry, because there may be two pairs sharing a site, and the chemistry of the other pair was NOT updated}
	 o1:=ich[l1]; o2:=ich[l2]; o3:=ich[k1]; o4:=ich[k2];
	 if found then found:=(chbin(pb1,o1,o2) and chbin(pb2,o3,o4));
	 if ((l mod 1000)=0) then writeln('No swappable pairs found after ',l:1,' attempts...');
      until found;
      nflp:=nflp+1; flptot:=flptot+1;
      {begin
      write('S2P:  ',l1:1,',',l2:1,' ',ich[l1]:1,' ',ich[l2]:1,' -> ');
      writeln(k1:1,',',k2:1,' ',ich[k1]:1,' ',ich[k2]:1,'  ');
      end;}
      {get u0}
      nppo:=4;
      ppo[1]:=l1; ppo[2]:=l2; ppo[3]:=k1; ppo[4]:=k2; {needed in site_ene_list}
      xp:=0.0;
      site_ene_list(l1,nppo,ppo,x,xp); u0:=x;
      site_ene_list(l2,nppo,ppo,x,xp); u0:=u0+x;
      site_ene_list(k1,nppo,ppo,x,xp); u0:=u0+x;
      site_ene_list(k2,nppo,ppo,x,xp); u0:=u0+x;
      {writeln(xp:1:5);}
      u0:=u0-xp/2.0;
      {save old chemistry}
      {another choice: k1->l1 k2->l2 or the other way?}
      x:=randu(seed);
      if (x>0.5) then begin
	 ich[l1]:=o3;
	 ich[l2]:=o4;
      end else begin
	 ich[l1]:=o4;
	 ich[l2]:=o3;
      end;
      x:=randu(seed);
      if (x>0.5) then begin
	 ich[k1]:=o1;
	 ich[k2]:=o2;
      end else begin
	 ich[k1]:=o2;
	 ich[k2]:=o1;
      end;
      if debug[5] then write(' ',o1:1,' ',o2:1,' ',o3:1,' ',o4:1,' -> ',ich[l1]:1,' ',ich[l2]:1,' ',ich[k1]:1,' ',ich[k2]:1,' ');
      {get u1}
      xp:=0.0;
      site_ene_list(l1,nppo,ppo,x,xp); u1:=x;
      site_ene_list(l2,nppo,ppo,x,xp); u1:=u1+x;
      site_ene_list(k1,nppo,ppo,x,xp); u1:=u1+x;
      site_ene_list(k2,nppo,ppo,x,xp); u1:=u1+x;
      {writeln(xp:1:5);}
      u1:=u1-xp/2.0;
      if debug[5] then write('dU=',(u1-u0):1:6,' ');
      x:=beta[kkl];
      du0:=x*(u1-u0);
      eavr[kkl,0]:=eavr[kkl,0]+1.0;
      eavr[kkl,lab]:=eavr[kkl,lab]+1.0;
      if (du0>eeps) then begin
	 x:=exp(-du0);
	 y:=randu(seed);
	 if (y>x) then begin {reject flip}
	    if debug[5] then write(' rejected');
	    ich[l1]:=o1; ich[l2]:=o2; ich[k1]:=o3; ich[k2]:=o4;
	    nflp:=nflp-1;
	    eavr[kkl,0]:=eavr[kkl,0]-1.0;
	    eavr[kkl,lab]:=eavr[kkl,lab]-1.0;
	 end else begin {accept}
	    if debug[5] then write('accepted1');
	    du:=du+u1-u0;
	    site_occ;
	    mystat(kkl,lab,du+epot);
	    pop1[p1]:=po2;
	    pop2[p2]:=po1;
	    if ((kk mod save_e)=0) then begin
	       write(recu,(u1-u0):1:6,' ');
	       write(recu,'   ',kkl:1,' ',lab:1,'   ');
	       writeln(recu,ich[l1]:1,' ',ich[l2]:1,' ',ich[k1]:1,' ',ich[k2]:1);
	    end;
	 end;
      end else begin {accept}
	 if debug[5] then write('accepted2');
	 du:=du+u1-u0;
	 site_occ;
	 mystat(kkl,lab,du+epot);
	 pop1[p1]:=po2;
	 pop2[p2]:=po1;
	 if ((kk mod save_e)=0) then begin
	    write(recu,(u1-u0):1:6,' ');
	    write(recu,'   ',kkl:1,' ',lab:1,'   ');
	    writeln(recu,ich[l1]:1,' ',ich[l2]:1,' ',ich[k1]:1,' ',ich[k2]:1);
	 end;
      end;
      x:=(epot+du)/natoc;
      if debug[5] then begin 
	 {for i:=1 to sites do n_a[pl[i,1],pl[i,2],1]:=ich[i];
	  UPPN(x); writeln('E-E_test = ',(epot+du-x):1:5);}
	 writeln(' E=',x:1:5);
      end;
      if (x<emin) then begin
	 list2atoms;
	 export_config_if(emin,x,beta[kkl],sname[lab]);
      end;
   end;
   list2atoms;
   if (flptot>0) then rate:=nflp/flptot else rate:=0.0;
   if debug[5] then writeln((sites-natoc):1,' vacancies.');
   epot:=epot+du;
   if debug[5] then begin
      write('Swap U rate ',(epot)/natoc:1:6,rate:8:4);
      UPPN(epot); writeln('  UPPN: ',epot/natoc:1:6);
   end;
end; { SWAP_two_pairs }
procedure SWAP_pair(kkl,mcs,lab,pb : integer;Var epot,rate:myreal);
{nn-tables, precalc energies...}
Var
   i,j,l,o,o1,o2     : integer;
   l1,l2,k1,k2	     : integer;
   kk,kk2	     : integer;
   ar,x,y	     : myreal;
   u0,u1,du0,du,ux   : myreal;
   v,r1,r2,w	     : v3;
   flptot,nflp	     : integer;
   swappable,identic : boolean;
   check	     : boolean;
   mfact	     : myreal;
   popa		     : array[1..nnpair] of integer; {pointers to useful pairs}
   npop		     : integer;
   found	     : boolean;
   ll		     : i10;
   
function chbin(p,ich1,ich2 : integer):boolean;
begin
   chbin:=false;
   if ((ich1=bondch[p,1]) and (ich2=bondch[p,2])) then chbin:=true;
   if ((ich1=bondch[p,2]) and (ich2=bondch[p,1])) then chbin:=true;
   {-1 is wildcard:}
   if ((ich1=bondch[p,1]) and (bondch[p,2]=-1)) then chbin:=true;
   if ((ich1=bondch[p,2]) and (bondch[p,1]=-1)) then chbin:=true;
   if ((ich2=bondch[p,1]) and (bondch[p,2]=-1)) then chbin:=true;
   if ((ich2=bondch[p,2]) and (bondch[p,1]=-1)) then chbin:=true;
   if ((bondch[p,1]=-1) and (bondch[p,2]=-1)) then chbin:=true;
end; { chbin }
function chbin_OR(p,ich1,ich2 : integer):boolean;
begin
   chbin_OR:=false;
   if (((ich1=bondch[p,1]) or (ich1=bondch[p,2])) AND
       ((ich2=bondch[p,1]) or (ich2=bondch[p,2]))) then chbin_OR:=true;
   {-1 is wildcard:}
   if (((ich1=bondch[p,1]) or (bondch[p,1]=-1)) AND
       ((ich2=bondch[p,2]) or (bondch[p,2]=-1))) then chbin_OR:=true;
   if (((ich2=bondch[p,1]) or (bondch[p,1]=-1)) AND
       ((ich1=bondch[p,2]) or (bondch[p,2]=-1))) then chbin_OR:=true;
   if ((bondch[p,1]=-1) and (bondch[p,2]=-1)) then chbin_OR:=true;
end; { chbin }
begin
   occupied_sites(natoc);
   if debug[5] then write('E_ini, Na  ',epot/natoc:1:5,' ',natoc:1);
   {list of bonds of type=pb}
   npop:=0; for i:=1 to npair do if bpair[i,pb] then begin
      npop:=npop+1; popa[npop]:=i;
   end;
   {in the destroy-recreate-list reshuffling method (as opposed to
    swap-site-list method) we have to demand than only the atoms belonging
    to the same ORBIT may swap; otherwise destruction-reconstruction of a motif
    bound to an object does not conserve chemistry:}
   check:=false; {will check equality of orbit# and layer# for a swap}
   {lets go for it}
   nflp:=0; flptot:=0;
   du:=0.0;
   kk2:=mcs*npop;
   if debug[5] then writeln('mcs npop nmax : ',mcs:1,' ',npop:1,' ',nmax:1);
   if (npop>1) then for kk:=1 to kk2 do begin
      eavnn[kkl,0]:=eavnn[kkl,0]+1.0;
      eavnn[kkl,lab]:=eavnn[kkl,lab]+1.0;
      x:=randu(seed);
      o1:=trunc(x*npop)+1;
      if (o1>npop) then o1:=npop;
      o2:=popa[o1];
      l1:=lpair[o2,1];
      l2:=lpair[o2,2];
      if ((ich[l1]<>ich[l2]) and (chbin(pb,ich[l1],ich[l2]))) then begin
	 nflp:=nflp+1; flptot:=flptot+1;
	 if (debug[5]) then write('SWAP_pair: ',l1:1,'/',l2:1,' ',ich[l1]:1,'/',ich[l2]:1,' -> ');
	 {get u0}
	 ux:=0.0;
	 ll[1]:=l2;
	 site_ene_list(l1,1,ll,x,ux); u0:=x;
	 ll[1]:=l1;
	 site_ene_list(l2,1,ll,x,ux); u0:=u0+x-ux/2.0;
	 {swap chemistries}
	 i:=ich[l1];
	 ich[l1]:=ich[l2];
	 ich[l2]:=i;
	 {get u1}
	 ux:=0.0;
	 ll[1]:=l2;
	 site_ene_list(l1,1,ll,x,ux); u1:=x;
	 ll[1]:=l1;
	 site_ene_list(l2,1,ll,x,ux); u1:=u1+x-ux/2.0;
	 if debug[5] then write('dU=',(u1-u0):1:6,' ');
	 x:=beta[kkl];
	 du0:=x*(u1-u0);
	 eavr[kkl,0]:=eavr[kkl,0]+1.0;
	 eavr[kkl,lab]:=eavr[kkl,lab]+1.0;
	 if (du0>eeps) then begin
	    x:=exp(-du0);
	    y:=randu(seed);
	    if (y>x) then begin {reject flip}
	       if debug[5] then write(' rejected');
	       i:=ich[l1];
	       ich[l1]:=ich[l2];
	       ich[l2]:=i;
	       nflp:=nflp-1;
	       eavr[kkl,0]:=eavr[kkl,0]-1.0;
	       eavr[kkl,lab]:=eavr[kkl,lab]-1.0;
	    end else begin {accept}
	       if debug[5] then write('accepted1');
	       du:=du+u1-u0;
	       site_occ;
	       mystat(kkl,lab,du+epot);
	       if ((kk mod save_e)=0) then begin
		  writeln(recu,(u1-u0):1:6,' ',kkl:1,' ',lab:1,'   ',ich[l1]:1,' ',ich[l2]:1);
	       end;
	    end;
	 end else begin {accept}
	    if debug[5] then write('accepted2');
	    du:=du+u1-u0;
	    site_occ;
	    mystat(kkl,lab,du+epot);
	    if ((kk mod save_e)=0) then begin
	       writeln(recu,(u1-u0):1:6,' ',kkl:1,' ',lab:1,'   ',ich[l1]:1,' ',ich[l2]:1);
	    end;
	 end;
	 x:=(epot+du);
	 if debug[5] then writeln(' E=',x/natoc:1:6);
	 if (x<emin) then begin
	    list2atoms;
	    export_config_if(emin,x,beta[kkl],sname[lab]);
	 end;
      end;
   end;
   list2atoms;
   if (flptot>0) then rate:=nflp/flptot else rate:=0.0;
   if debug[5] then writeln((sites-natoc):1,' vacancies.');
   epot:=epot+du;
   if (debug[5] or debug[3]) then begin
      write('Swap U rate ',(epot)/natoc:1:6,rate:8:4);
      UPPN(epot); writeln('  UPPN: ',epot/natoc:1:6);
   end;
end; { swap_pair }
procedure AtomFlip(kkl,mcs,lab:integer;Var epot,rate:myreal);
{nn-tables, precalc energies...}
Var
   i,j,l,o,nppo	: integer;
   p1,l1,l2,kk	: integer;
   ar,x,y,xp	: myreal;
   u0,u1,du0,du	: myreal;
   v,r1,r2,w	: v3;
   flptot,nflp	: integer;
   poo		: array[1..nsites,1..2] of integer;
   npo,mu	: integer;
   ch10,ch10x	: i10;
   s10		: array[1..10] of i10;
   flippable	: boolean;
   ppo		: i10;
   rnm,rn5	: integer;
begin
   for i:=1 to 10 do for j:=1 to 10 do s10[i,j]:=0;
   {first, allocate the orbits subject to flipping}
   npo:=0; mu:=0;
   for i:=1 to sites do begin
      mu:=mu+1;
      if (mu=1) then begin
	 npo:=npo+1;
	 poo[npo,1]:=i;
	 poo[npo,2]:=omu[och[i]];
      end;
      if (mu=omu[och[i]]) then mu:=0;
   end;
   if debug[6] then writeln('AtomFlip: ',sites:1,' sites, max. ',nmax:1,' neighbors, ',npo:1,' node-orbits.');
   occupied_sites(natoc);
   {lets go for it}
   nflp:=0; flptot:=0;
   du:=0.0;
   for kk:=1 to mcs*npo do begin
      eavnn[kkl,0]:=eavnn[kkl,0]+1.0;
      eavnn[kkl,lab]:=eavnn[kkl,lab]+1.0;
      x:=randu(seed);
      o:=trunc(x*npo)+1;
      if o>npo then o:=npo;
      l1:=poo[o,1]; {pointer}
      l2:=poo[o,2]; {multiplicity}
      flippable:=false;
      flptot:=flptot+1;
      {for i:=0 to l2-1 do write(ich[l1+i]:1,' '); writeln;}
      for i:=1 to l2-1 do if (ich[l1+i]<>ich[l1+i-1]) then flippable:=true;
      if flippable then begin
	 nflp:=nflp+1;
	 {get u0}
	 {site_ene_list is needed to correct for double-counted bonds}
	 nppo:=0; for i:=0 to l2-1 do if (ich[l1+i]>0) then begin
	    nppo:=nppo+1; ppo[nppo]:=l1+i;
	 end;
	 u0:=0.0; xp:=0.0; for i:=1 to nppo do begin
	    site_ene_list(ppo[i],nppo,ppo,x,xp); u0:=u0+x;
	 end; u0:=u0-xp/2;
	 {u0:=0.0; for i:=0 to l2-1 do if (ich[l1+i]>0) then begin
	 site_ene(l1+i,x); u0:=u0+x;
	 end;}
	 if debug[6] then write(u0:1:6,' ');
	 {swap: note we prevent p1=1 (self-mapping)}
	 if (l2=2) then p1:=2 else begin
	    x:=randu(seed)*5; {random number for 5-fold rotation [1-5]}
	    rn5:=round(x+0.5);
	    if (rn5=0) then rn5:=1; if (rn5>5) then rn5:=5;
	    x:=randu(seed)*2; {random number for mirror: [1-2]}
	    rnm:=round(x+0.5);
	    if (rnm=0) then rnm:=1; if (rnm>2) then rnm:=2;
	    {x:=randu(seed); p1:=trunc(x*(l2-1))+2;
	    if (p1>l2) then p1:=l2;} {shouldnt happen}
	 end;
	 {s10[l2,p1]:=s10[l2,p1]+1;}
	 for o:=0 to l2-1 do ch10[o+1]:=ich[l1+o];
	 ch10x:=ch10; SHUFFLE(l2,rn5,rnm,ch10x);
	 {for o:=1 to l2 do ch10x[gsym[l2,1,o]]:=ch10[gsym[l2,p1,o]];}
	 for o:=0 to l2-1 do ich[l1+o]:=ch10x[o+1];
	 if debug[6] then begin
	    writeln(l1:1,' ',l2:1,' ',och[l1]:1,' ',p1:1);
	    for o:=1 to l2 do write(ch10[o]:3); writeln;
	    for o:=1 to l2 do write(ch10x[o]:3); writeln;
	 end;
	 {get u1}
	 nppo:=0; for i:=0 to l2-1 do if (ich[l1+i]>0) then begin
	    nppo:=nppo+1; ppo[nppo]:=l1+i;
	 end;
	 u1:=0.0; xp:=0.0; for i:=1 to nppo do begin
	    site_ene_list(ppo[i],nppo,ppo,x,xp); u1:=u1+x;
	 end; u1:=u1-xp/2;
	 {u1:=0.0; for i:=0 to l2-1 do if (ich[l1+i]>0) then begin
	 site_ene(l1+i,x); u1:=u1+x;
	 end;}
	 if debug[6] then write(u1:1:6,' ');
	 x:=beta[kkl];
	 du0:=x*(u1-u0);
	 eavr[kkl,0]:=eavr[kkl,0]+1.0;
	 eavr[kkl,lab]:=eavr[kkl,lab]+1.0;
	 if (du0>eeps) then begin
	    x:=exp(-du0);
	    y:=randu(seed);
	    if y>x then begin {reject flip}
	       if debug[6] then writeln(' rejected');
	       for o:=0 to l2-1 do ich[l1+o]:=ch10[o+1];
	       nflp:=nflp-1;
	       eavr[kkl,0]:=eavr[kkl,0]-1.0;
	       eavr[kkl,lab]:=eavr[kkl,lab]-1.0;
	    end else begin {accept}
	       if debug[6] then writeln('accepted1');
	       du:=du+u1-u0;
	       site_occ;
	       mystat(kkl,lab,du+epot);
	       if ((kk mod save_e)=0) then begin
		  writeln(recu,du:1:6,' ',kkl:1,' ',lab:1);
	       end;
	    end;
	 end else begin {accept}
	    if debug[6] then writeln('accepted2');
	    du:=du+u1-u0;
	    site_occ;
	    mystat(kkl,lab,du+epot);
	    if ((kk mod save_e)=0) then begin
	       writeln(recu,du:1:6,' ',kkl:1,' ',lab:1);
	    end;
	 end;
	 x:=(epot+du)/natoc;
	 if (x<emin) then begin
	    list2atoms;
	    export_config_if(emin,x,beta[kkl],sname[lab]);
	 end;
      end;
   end;
   list2atoms;
   rate:=nflp/flptot;
   if debug[6] then writeln((sites-natoc):1,' vacancies.');
   epot:=epot+du;
   if debug[6] then begin
      write('AtomFlip U rate ',(epot)/natoc:1:6,rate:8:4);
      UPPN(epot); writeln('  UPPN: ',epot/natoc:1:6);
   end;
   {for i:=1 to 10 do begin
    k:=0; for o:=1 to 10 do k:=k+s10[i,o];
    if k>0 then begin
    write(i:2); for o:=1 to i do write(s10[i,o]/nflp:8:4); writeln;
      end;
   end;}
end; { AtomFlip }

procedure random_init;
var
   i,l,kk,l1,l2  : integer;
   sitesx,iat,lbl : integer;
   x,rate	  : myreal;
begin
   occupied_sites(natoc);
   {Stage one: remove some atoms when n_occ>n_nom}
   for o:=1 to nnele do if (mele_oc[o]>mele_nom[o]) then begin
      l1:=mele_oc[o]-mele_nom[o];
      repeat
	 x:=randu(seed);
	 i:=trunc(x*sites)+1; if i>sites then i:=sites;
	 if (ich[i]=ele[o]) then begin
	    ich[i]:=0; l1:=l1-1;
	    mele_oc[o]:=mele_oc[o]-1; natoc:=natoc-1;
	 end;
      until (l1=0);
   end;
   if debug[4] then write('Reinitializing chemistry: ');
   {input consistency check}
   iat:=0; for o:=1 to nnele do iat:=iat+mele_nom[o];
   if iat>sites then begin
      writeln('Inconsistent input: atom# exceeds site#: ',iat:1,'>',sites:1);
      halt;
   end;
   sitesx:=0;
   for i:=1 to sites do  {vacancy and count frozen sites}
      if not(s_frozen[i]) then begin nfree:=nfree+1
      end else sitesx:=sitesx+1;
   for o:=1 to nnele do begin {counts of frozen sites}
      mele_fro[o]:=0;
      for i:=1 to sites do
	 if (s_frozen[i] and (ich[i]=ele[o])) then mele_fro[o]:=mele_fro[o]+1;
   end;
   j:=0; for o:=1 to nnele do j:=j+mele_fro[o];
   if (j<>sitesx) then begin
      writeln('Inconsistent chemistry reinit (frozen orbits): ');
      for o:=1 to nnele do write(ele[o]:1,' ',mele_fro[o]:1,'  '); writeln;
      writeln('Sum of the above not equal ',sitesx:1);
      halt;
   end;
   for o:=1 to nnele do if (mele_fro[o]>mele_nom[o]) then begin
      writeln('d04.inp : set N_[',ele[o]:1,']=',mele_fro[o]:1,' or x_[',ele[o]:1,']>=',mele_fro[o]/natnom:1:4,' (frozen-orbit constraint).');
      halt;
   end;
   if debug[4] then begin
      write('Frozen sites: ');
      for o:=1 to nnele do write(ele[o]:1,' ',mele_fro[o]:1,'   '); writeln;
   end;
   {writeln('Free sites: ',nfree:1);}
   nchs:=0; for o:=1 to nnele do mele_buf[o]:=0;
   for o:=1 to nnele do if (mele_oc[o]<mele_nom[o]) then begin
      j:=0; repeat j:=j+1;
	 nchs:=nchs+1; chlist[nchs]:=ele[o];
	 mele_buf[o]:=mele_buf[o]+1;
	 {natoc:=natoc-1;}
      until ((j+mele_oc[o])=mele_nom[o]);
   end;
   if debug[4] then begin
      write('Buffered sites: ');
      for o:=1 to nnele do write(ele[o]:1,' ',mele_buf[o]:1,'   '); writeln;
      write('Occupied sites: ',natoc:1,'   ');
      for o:=1 to nnele do write(ele[o]:1,' ',mele_oc[o]:1,'   '); writeln;
   end;
   {check chlist:}
   for o:=1 to nnele do begin
      j:=0; for i:=1 to nchs do if ele[o]=chlist[i] then j:=j+1;
      if (mele_nom[o]-mele_oc[o]<>j) then begin
	 writeln('Error creating chlist:',ele[o]:1,' -> ',mele_oc[o]:1,'+',j:1,'<>',mele_nom[o]:1); halt;
      end;
   end;
   if debug[4] then writeln('Random-init (stage 2): adding ',(natlim-natoc):1,' atoms -> ');
   if (natlim>natoc) then begin
      {if ((nbond>0) and (mcs[loops+1]>0)) then enn_table;}
      epot:=0.0; rate:=0.0; {doesnt matter here}
      repeat
	 x:=randu(seed);
	 i:=trunc(x*nchs)+1; if i>nchs then i:=nchs;
	 repeat
	    x:=randu(seed);
	    l:=trunc(x*sites)+1; if l>sites then l:=sites;
	 until (not(s_frozen[l]) and (ich[l]=0));
	 ich[l]:=chlist[i];
	 l1:=pele[ich[l]];
	 natoc:=natoc+1;
	 {if mcs[loops+1]>0 then begin
	  for j:=1 to nbond do begin
	  lbl:=2; kk:=loops+1;
	  swap_pair(kk,mcs[kk],lbl,j,epot,rate);
	    end;
	 end;}
	 for j:=i to nchs-1 do chlist[j]:=chlist[j+1];
	 nchs:=nchs-1;
	 mele_oc[l1]:=mele_oc[l1]+1;
	 mele_buf[l1]:=mele_buf[l1]-1;
	 {writeln(l:1,' ',l1:1,' ',ich[l]:1,' ',mele_oc[l1]:1);}
      until ((nchs=0) or (natoc=natlim));
   end;
   if debug[4] then begin
      write('Occupied sites: ',natoc:1,'  ');
      for o:=1 to nnele do write(ele[o]:1,' ',mele_oc[o]:1,'   '); writeln;
   end;
   for i:=1 to sites do begin
      n_a[pl[i,1],pl[i,2],1]:=ich[i];
      add_chem_para(n_para[pl[i,1],pl[i,2]],ich[i]);
   end;
end; { random_init }


procedure LatticeGas(reinit: boolean;kkl,mcs,lab:integer;Var epot,rate:myreal);
{nn-tables, precalc energies...}
Var
   i,j,l,o	       : integer;
   i1,i2,l1,l2,kk,nflp : integer;
   ar,x,y	       : myreal;
   u0,u1,du0,du	       : myreal;
   u,ud,udd,chx	       : myreal;
   v,r1,r2,w	       : v3;
   flptot	       : integer;
   attempt_add	       : boolean;
   attempt_remove      : boolean;

begin
   nflp:=0; flptot:=0;
   du:=0.0;
   occupied_sites(natoc);
   for kk:=1 to mcs*sites do begin
      for i:=1 to nchs do if chlist[i]=0 then begin
	 writeln('chlist[',i:1,']=0 .'); halt;
      end;
      nflp:=nflp+1; flptot:=flptot+1;
      eavnn[kkl,0]:=eavnn[kkl,0]+1.0;
      eavnn[kkl,lab]:=eavnn[kkl,lab]+1.0;
      x:=randu(seed);
      l:=trunc(x*sites)+1; if l>sites then l:=sites;
      if not(s_frozen[l]) then begin
	 if debug[4] then write(l:1,' ',ich[l]:1,' : ');
	 attempt_remove:=(ich[l]>0);
	 attempt_add:=(ich[l]=0);
	 if attempt_remove then begin {remove atom?}
	    chx:=chpt[ich[l]];
	    site_ene(l,u0); u0:=-u0-u_self[ich[l]]-chx;
	    if debug[4] then write(' ch>0 ',u0:1:4,' ',chx:1:2);
	    du0:=beta[kkl]*u0;
	    if (du0>eeps) then begin
	       x:=exp(-du0);   
	       y:=randu(seed); if (y<eeps) then writeln('rnd ',y:1:7,' ',x:1:7);
	       if (y>x) then begin {reject flip}
		  nflp:=nflp-1;
		  if debug[4] then writeln(' ',nchs:1,' ',natoc:1,' rejected');
	       end else begin {accept u-incr.}
		  l1:=pele[ich[l]];
		  natoc:=natoc-1;
		  nchs:=nchs+1;
		  chlist[nchs]:=ich[l];
		  ich[l]:=0;
		  mele_fro[l1]:=mele_fro[l1]+1;
		  mele_buf[l1]:=mele_buf[l1]+1;
		  du:=du+u0;
		  if not(reinit) then begin
		     site_occ;
		     mystat(kkl,lab,du+epot);
		  end;
		  if debug[4] then writeln(' ',nchs:1,' ',natoc:1,' +accepted');
	       end;
	    end else begin {accept}
	       l1:=pele[ich[l]];
	       natoc:=natoc-1;
	       nchs:=nchs+1;
	       chlist[nchs]:=ich[l];
	       ich[l]:=0;
	       mele_fro[l1]:=mele_fro[l1]+1;
	       mele_buf[l1]:=mele_buf[l1]+1;
	       du:=du+u0;
	       if not(reinit) then begin
		  site_occ;
		  mystat(kkl,lab,du+epot);
	       end;
	       if debug[4] then writeln(' ',nchs:1,' ',natoc:1,' accepted');
	    end;
	 end;
	 if attempt_add and (nchs>0) then begin {add atom?}
	    if chlist[nchs]=0 then begin
	       writeln('last element in chlist = 0 . nchs=',nchs:1,' .');
	       halt;
	    end;
	    x:=randu(seed);
	    l1:=trunc(x*nchs)+1;
	    if (l1>nchs) then l1:=nchs;
	    if (l1<1) then l1:=1;
	    if chlist[l1]=0 then begin
	       writeln('Error : chlist[',l1:1,']=0 .'); halt;
	    end;
	    ich[l]:=chlist[l1];
	    chx:=chpt[ich[l]];
	    site_ene(l,u0); u0:=u0+u_self[ich[l]]+chx;
	    if debug[4] then write(' ',chlist[nchs]:1,' ',l1:1,' ',chlist[l1]:1,' ',u0:1:4,' ',chx:1:2);
	    du0:=beta[kkl]*u0;
	    if (du0>eeps) then begin
	       x:=exp(-du0);   
	       y:=randu(seed);
	       if (y>x) then begin {reject flip}
		  ich[l]:=0;
		  nflp:=nflp-1;
		  if debug[4] then writeln(' ',nchs:1,' ',natoc:1,' rejected');
	       end else begin {accept incr.}
		  if (l1<nchs) then
		     for i:=l1 to nchs-1 do chlist[i]:=chlist[i+1];
		  chlist[nchs]:=0;
		  nchs:=nchs-1;
		  natoc:=natoc+1;
		  l1:=pele[ich[l]];
		  mele_fro[l1]:=mele_fro[l1]-1;
		  mele_buf[l1]:=mele_buf[l1]-1;
		  du:=du+u0;
		  if not(reinit) then begin
		     site_occ;
		     mystat(kkl,lab,du+epot);
		  end;
		  if debug[4] then writeln(' ',nchs:1,' ',natoc:1,' +accepted');
	       end;
	    end else begin {accept}
	       for i:=l1 to nchs-1 do chlist[i]:=chlist[i+1];
	       chlist[nchs]:=0;
	       nchs:=nchs-1;
	       natoc:=natoc+1;
	       l1:=pele[ich[l]];
	       mele_fro[l1]:=mele_fro[l1]-1;
	       mele_buf[l1]:=mele_buf[l1]-1;
	       du:=du+u0;
	       if not(reinit) then begin
		  site_occ;
		  mystat(kkl,lab,du+epot);
	       end;
	       if debug[4] then writeln(' ',nchs:1,' ',natoc:1,' accepted');
	    end;
	 end;
      end;
   end;
   list2atoms;
   if reinit then begin
      occupied_sites(natoc);
      UPPN(epot);
   end else epot:=epot+du;
   rate:=nflp/flptot;
   if debug[4] then writeln((sites-natoc):1,' vacancies.');
end; { LatticeGas }
procedure init_averages;
var i,o,j,k : integer;
begin
   for i:=1 to tloops do
      for o:=0 to tact do begin {running averages of E}
	 eav[i,o]:=0.0; eavn[i,o]:=0.0; eavr[i,o]:=0.0; eav2[i,o]:=0.0;
	 eavnn[i,o]:=0.0;
      end;
      for j:=1 to t_orb do begin
	 o_al[i,j]:=0.0;
	 for k:=0 to z_ele do o_ch[i,j,k]:=0.0;
      end;
end; { init_averages }       
procedure d04_main;
Var
   i,j,k,l,o,p1,p2	 : integer;
   ii,jj,kk,ll,mm,dep	 : integer;
   imcs,lbl,ntw,nsw	 : integer;
   u,u0,u1		 : myreal;
   epot,e0,e1,x,y	 : myreal;
   rate,dr2		 : myreal;
   txf			 : text;
   elg,hcap		 : myreal;
   str3			 : string[3];
   reinit		 : boolean;   
   chin,rwf,rf,ueff,rmsd : myreal;
   chi3,echi3		 : v3;
   nsamples		 : integer;
   ss			 : mystring;
begin
   kboltz:=1.0/11604.5;
   if(nflipnu=0) then begin
      for i:=-tfla to tfla do
	 for j:=-tfla to tfla do flip_b[i,j]:=true;
   end else begin
      for i:=-tfla to tfla do
	 for j:=-tfla to tfla do flip_b[i,j]:=false;
   end;
   
   {flip_nu[nflipnu,1]:=dperp_offset(nu_offs,j);}
   {flip_nu[nflipnu,2]:=dperp_offset(nu_offs,k);}
   for i:=1 to nflipnu do begin
      {k:=dperp_offset(nu_offs,flip_nu[i,1]);}
      {l:=dperp_offset(nu_offs,flip_nu[i,2]);}
      k:=flip_nu[i,1];
      l:=flip_nu[i,2];
      flip_b[k,l]:=true;
      flip_b[l,k]:=true;
      if debug[2] then writeln('DPERP: flip ',k:1,' ',l:1,' allowed: ',flip_b[k,l]);
   end;
   if not(flavours) then begin
      for i:=-tfla to tfla do
	 for j:=-tfla to tfla do flip_b[i,j]:=true;
   end;
   ntocc:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do
	 for k:=1 to t_ele do n_oc[i,j,k]:=0;
   rewrite(recu,'u.rep');
   if deco then rewrite(txta,'copas.out');
   eeps:=0.0; {dE in U/kT units that makes any difference}
   epot:=0.0;
   elg:=0.0; rate:=0.0;
   reinit:=false;
   if(hiter<hiterw)then hiterw:=hiter;
   if(hiter>=hiterw)then nsamples:=loops else nsamples:=loops*hiterw;
   if deco then begin
      occupied_sites(natt); {writeln(natt:1,' atoms.');}
      for i:=1 to nnele do begin
	 {writeln('DBG_CHVAR_BEG ',export_err,' ',i:1,' ',ele[i]:1,' ',mele_oc[i]:1);}
	 chnom[ele[i]]:=chcur[ele[i]];
      end;
      if save_xyzs then begin
	 writeln('SAVE_XYZS is ON: export ',hiter*viter*loops:1,' samples (after each horizontal cycle) -> xyzs.d04.');
	 rewrite(txyzs,'xyzs.d04');
	 writeln(txyzs,viter*hiter:1,' ',loops:1,' ',natt:1,'   # viter*hiter, T-loops, NAT ');
	 for i:=1 to 3 do begin
	    for o:=1 to 3 do write(txyzs,bl3[i,o]:10:5); writeln(txyzs);
	 end;
      end;
      if (deco) then begin
	 if(flip_man) then
	    write(txta,tnodes:1,' ',multi:1,' ',(nseq+1):1) else
	       write(txta,tnodes:1,' ',multi:1,' ',nsamples*viter:1);
	 if(instances) then writeln(txta,'  # read_tvl=1') else writeln(txta);
	 for i:=1 to 2 do begin
	    for o:=1 to 5 do write(txta,bai[i,o]:5); writeln(txta);
	 end;
      end;
      if (cheminp) then begin
	 atom_list;
	 random_init;
	 for o:=1 to nnele do chnom[ele[o]]:=chcur[ele[o]];
	 if debug[11] then write('Random-Init: ');
	 occupied_sites(natt); write(natt:1,' atoms, ');
	 UPPN(epot);
      end else UPPN(epot);
      if (chvar) then begin
	 chempot(cmp);
	 epot:=epot+cmp;
      end;
      writeln('E_ini=',epot/natt:1:6,' [eV/atom]');
      if debug[12] then export_config('copa.ini');
   end;
   if not(deco) then begin
      init_th(epot);
      natt:=tnodes;
      writeln('E_ini=',epot/tnodes:1:6);
   end;
   {export tiling after finishing a loop}
   if flip_replay then begin
      rewrite(txti,'tilings.mon');
      write(txti,tnodes:1,' ',multi:1,' ',(til2-til1+1):1,'    ');
      for o:=til1 to til2 do write(txti,fene[o]:1:3,' ');
   end else begin
      rewrite(txti,'tilings.out');
      write(txti,tnodes:1,' ',multi:1,' ',nsamples*hiter*viter:1,'    ');
      {for o:=1 to loops do write(txti,beta[o]:1:2,' ');}
   end;
   writeln(txti,'  # nodes loops beta/dU');
   for i:=1 to 2 do begin
      for o:=1 to 5 do write(txti,bai[i,o]:5); writeln(txti);
   end;
   if (flip_rec) then begin
      rewrite(txfout,'flip.rec');
      for i:=1 to 2 do begin
	 for o:=1 to 5 do write(txfout,bai[i,o]:5); writeln(txfout);
      end;
      writeln(txfout,tnodes:1,' ',multi:1,' ',deconame,' flip_record');
      export_tiling_body(txfout);
      writeln(txfout,couple:1,'  couple');
   end;
   init_averages;
   if table_reinit then enn_table;
   {if ((relmcs>0) or (dodfit)) then nntable(rcplus);}
   writeln(recu,epot:1:6,' ',tnodes:1,' ',natt:1,'  # E_ini NODES NAT');
   for ivi:=1 to viter do begin {vertical iteration}
      for kk:=1 to loops do begin {T-cycle}
	 nflp:=0; ntflp:=0; 
	 for o:=-tfla to tfla do begin
	    nflx[o]:=0;
	 end;
	 for o:=1 to t_obj do noflx[o]:=0;
	 for ihi:=1 to hiter do begin {horizontal iteration}
	    ntw:=0; nsw:=0;
	    if (hiter-ihi<hiterw) then write_stout:=true else
	       write_stout:=false;
	    for ll:=1 to nschemes do begin
	       if (write_stout and not(flip_man or flip_replay)) then write("loop ",ivi:1,":",kk:1,":",ihi:1,":",ll:1," ");
	       {RELAXATION}
	       if ((sch[ll,1]=7) and (relmcs>0)) then begin
		  nntable(rcplus);
		  och_status;
		  RELAX_PARA(sch[ll,2],cueff,xgamma,dramax,dr2,epot); {7}
		  if(cueff>0.00001) then begin
		     ueff:=epot/natt;
		     UPPAV(epot,se);
		  end;
		  if write_stout then begin
		     write('   ',epot/natt:1:5,' ',sname[7]);
		     if(cueff>0.00001) then write(' ueff=',ueff:1:5);
		     writeln;
		  end;
		  if table_reinit then enn_table;
	       end;
	       {RELAXATION}
	       if ((sch[ll,1]=9) and (sch[ll,2]>0)) then begin
		  nntable(rcplus);
		  CMC(sch[ll,2]*mcs[kk],kk,cueff,sigmc,rate,epot,ueff,rmsd); {9}
		  if write_stout then
		     writeln(beta[kk]:1:1,' ',epot/natt:1:5,' ',rate:1:5,' ',sname[9],': ueff=',ueff:1:4,' rmsd=',rmsd:1:3);
		  if table_reinit then enn_table;
	       end;
	       {DFIT}
	       if (sch[ll,1]=8) then begin
		  nntable(rcplus);
		  for mm:=1 to nfita do begin
		     if ((fita[mm,1]<0) or (fita[mm,2]<0)) then begin
			if (fita[mm,1]<0) then randomize_params(1); {positions}
			if (fita[mm,2]<0) then randomize_params(0); {dw factors}
		     end else begin
			lambdafit:=lambdas[kk];
			etarg:=etargs[kk];
			DFIT(fita[mm,1],fita[mm,2],ueff,rf,rwf,chi3);
			if(fita[mm,1]>1) then echi3:=chi3;
			chin:=chi3[1];
		     end;
		  end;
		  chi3[3]:=echi3[3];
		  chin:=chi3[1];
		  if(cueff>1e-10) then UPPAV(epot,se) else epot:=ueff;
		  writeln(kk:1,' ',ihi:1,' ',etarg:1:5,' ',ueff/natoms:1:5,' ',epot/natoms:1:5,' ',chin:1:3,' ',rwf:1:4,' ',rf:1:4,' #Et Ueff E CH2 RWF RF');
		  if(chin<chimin)then begin
		     ss:='Exporter: chisq-min ';
		     paradump(txpardf,kk,ihi,epot,ueff,rf,rwf,chi3);
		     chimin:=chin;
		  end;
		  write(txrepdf,kk:1,' ',ihi:1,' ',ueff/natoms:1:5,' ',epot/natoms:1:5,' ',rwf:1:4,' ',rf:1:4);
		  for o:=1 to 3 do write(txrepdf,' ',chi3[o]:1:3); writeln(txrepdf);
		  
		  {if write_stout then begin
		   write('   ',epot/natt:1:5,' ',sname[8]);
		   writeln(': chi=',chin:1:2,' rwf=',rwf:1:4,' rf=',rf:1:4,' mfit=',mfit:1,' mfitp=',mfitp:1);
		  end;}
		  if table_reinit then enn_table;
	       end;
	       {LATTICE GAS}
	       if ((sch[ll,1]=1) and (lmcs>0) and (lgchp)) then begin
		  lbl:=1; 
		  LATTICEGAS(reinit,kk,sch[ll,2]*mcs[kk],lbl,epot,rate); {1}
		  if write_stout then begin
		     write(beta[kk]:1:1,' ',epot/natt:1:5,' ',rate:1:5,' ',sname[lbl],'  nchs=',nchs:1,' ');
		     for o:=1 to nnele do begin
			write(ele[o]:1,' ',mele_buf[o]:1,'  ');
		     end; writeln;
		  end;
	       end;
	       {SWAP}
	       if (sch[ll,1]=2) then begin
		  nsw:=nsw+1;
		  if (npp[nsw]>0) then begin
		     nflp:=0; ntflp:=0; lbl:=2;
		     SWAP_PAIR(kk,sch[ll,2]*mcs[kk],lbl,nsw,epot,rate); {2}
		     occupied_sites(natt);
		     if write_stout then writeln(beta[kk]:1:1,' ',epot/natt:1:5,' ',rate:1:5,' ',sname[lbl],' (',nsw:1,') Nb=',npp[nsw]:1,' Na=',natt:1);
		     if debug[4] then enn_table;
		  end else if write_stout then writeln(sname[lbl],'(',nsw:1,'): 0 pairs.');
	       end;
	       {SWAP2}
	       if ((sch[ll,1]=6) and (smcs2>0) and (nbond>0)) then begin
		  ntw:=ntw+1;
		  nflp:=0; ntflp:=0; lbl:=6;
		  SWAP_TWO_PAIRS(kk,sch[ll,2]*mcs[kk],lbl,twopa[ntw,1],twopa[ntw,2],epot,rate); {5}
		  if write_stout then writeln(beta[kk]:1:1,' ',epot/natt:8:5,rate:1:5,' ',sname[lbl],'(',ntw:1,') ',twopa[ntw,1]:1,' ',twopa[ntw,2]:1,' Na=',natt:1);
	       end;
	       {ORBIT-FLIP (GFLIP)}
	       if (sch[ll,1]=3) then begin
		  nflp:=0; ntflp:=0; lbl:=3;
		  ATOMFLIP(kk,sch[ll,2]*mcs[kk],lbl,epot,rate); {3}
		  occupied_sites(natt);
		  if write_stout then writeln(beta[kk]:1:1,' ',epot/natt:8:5,' ',rate:1:5,' ',sname[lbl],' Na=',natt:1);
	       end;
	       {OBJECT-FLIP}
	       if (sch[ll,1]=4) then begin
		  ntflp:=0; nflp:=0; lbl:=4;
		  if not(deco) then begin natt:=tnodes; end;
		  fobj:=1; savesigma:=sigmaflip; sigmaflip:=false;
		  MC_TILING(fobj,couple,vpin,level,tlevel,kk,sch[ll,2]*mcs[kk],lbl,epot,rate); {4}
		  sigmaflip:=savesigma;
		  if (write_stout and not(flip_man or flip_replay)) then
		  begin
		     write(beta[kk]:1:1,' ',epot/natt:1:5,' ',rate:1:5,' ',sname[lbl],':');
		     for o:=1 to n_obj do if (noflx[o]>0) then write(' ',obj_name[o],'(',noflx[o]/ntflp:1:4,')');
		     for o:=1 to n_obj do noflx[o]:=0;
		     if deco then begin
			occupied_sites(natt);
			write(' Na=',natt:1);
		     end;
		     if not(deco) then begin
			for o:=1 to n_objch do write(' ',obj_sta[o]:1);
		     end;
		     writeln;
		  end;
		  if table_reinit then enn_table;
		  if ((relmcs>0) or (dodfit)) then nntable(rcplus);
	       end;
	       {INSTANCE-FLIP}
	       if (sch[ll,1]=10) then begin
		  ntflp:=0; nflp:=0; lbl:=10;
		  if not(deco) then begin natt:=tnodes; end;
		  for o:=-ttvl to 0 do ntvl[o]:=0;
		  for o:=1 to tnodes do ntvl[tvl[o]]:=ntvl[tvl[o]]+1;
		  fobj:=2;
		  MC_TILING(fobj,couple,vpin,level,tlevel,kk,sch[ll,2]*mcs[kk],lbl,epot,rate); {10}
		  for o:=-ttvl to 0 do ntvl[o]:=0;
		  for o:=1 to tnodes do ntvl[tvl[o]]:=ntvl[tvl[o]]+1;
		  if (write_stout and not(flip_man or flip_replay)) then
		  begin
		     write(beta[kk]:1:1,' ',epot/natt:1:5,' ',rate:1:5,' ',sname[lbl],':');
		     for o:=-ttvl+1 to -1 do if (ntvl[o]>0) then write(' inst',o:1,'(',ntvl[o]:1,')');
		     if deco then begin
			occupied_sites(natt);
			write(' Na=',natt:1);
		     end;
		     if not(deco) then begin
			for o:=1 to n_objch do write(' ',obj_sta[o]:1);
		     end;
		     writeln;
		  end;
		  if table_reinit then enn_table;
		  if ((relmcs>0) or (dodfit)) then nntable(rcplus);
	       end;
	       {NODE-FLIP (TILING)}
	       if (sch[ll,1]=5) then begin
		  nflp:=0; ntflp:=0; lbl:=5;
		  if not(deco) then begin natt:=tnodes; end;
		  fobj:=0;
		  {writeln('sigmaflip ',sigmaflip);}
		  MC_TILING(fobj,couple,vpin,level,tlevel,kk,sch[ll,2]*mcs[kk],lbl,epot,rate); {4}
		  if (write_stout and not(flip_man or flip_replay)) then
		  begin
		     {symmetrize running info with respect to flavours}
		     for o:=-tfla to -1 do if (nflx[o]>0) then begin
			nflx[-o]:=nflx[-o]+nflx[o];
			nflx[o]:=0;
		     end;
		     write(beta[kk]:1:1,' ',epot/natt:1:5,' ',rate:1:5,' ',sname[lbl]);
		     for o:=-tfla to tfla do if (nflx[o]>0) then write(' ',o:1,'(',nflx[o]/ntflp:1:4,')');
		     for o:=-tfla to tfla do nflx[o]:=0;
		     if deco then begin
			occupied_sites(natt);
			write(' Na=',natt:1);
		     end;
		     if not(deco) then begin
			for o:=1 to n_objch do write(' ',obj_sta[o]:1);
		     end;
		     writeln;
		  end;
		  if table_reinit then enn_table;
		  if ((relmcs>0) or (dodfit)) then nntable(rcplus);
	       end;
	       occupied_sites(natt);
	       if save_xyzs then begin
		  writeln(txyzs,beta[kk]:1:3,'  beta');
		  export_xyz_body(txyzs);
	       end;
	    end; {types of MC}
	    if (not(flip_debug or flip_replay) and write_stout) then begin
	       export_tiling_body(txti);
	       if (deco) then begin
		  {writeln('saving ',ivi:1);}
		  writeln(txta,beta[kk]:1:1,' ',epot/natt:1:6);
		  export_tiling_body(txta);
		  export_params(txta);
		  export_tiling_atoms(txta);
	       end;
	    end;
	 end; {horizontal iteration}
      end; {T-cycle}
   end; {vertical iteration}
   if not(deco) then natt:=tnodes;
   rewrite(txf,'hc.out');
   writeln(txf,'UNITS: beta [1/eV], heat capacity HC [10^-6*ev/K], U [ev/at]');
   write(txf,' beta    U/atom     TOTAL');
   for o:=1 to 5 do if (eavn[1,o]>0) then write(txf,sname[o]:11); writeln(txf);
   for kk:=1 to loops do begin
      if (eavn[kk,0]>0) then write(txf,beta[kk]:5:1,eav[kk,0]/eavn[kk,0]/natt:10:5);
      for o:=0 to tact do if (eavn[kk,o]>0) then begin
	 hcap:=beta[kk]*beta[kk]*(eav2[kk,o]/eavn[kk,o]-eav[kk,o]*eav[kk,o]/eavn[kk,o]/eavn[kk,o])/natt*kboltz*1000000.0;
	 write(txf,hcap:11:3);
      end;
      if (eavn[kk,0]>0) then writeln(txf);
   end;
   close(txf);
   if (test_glob=1) then begin
      cmp:=0;
      if deco then begin
	 UPPN(e1);
	 if (chvar) then chempot(cmp);
	 e0:=epot/natt;
	 e1:=(e1+cmp)/natt;
      end;
      if not(deco) then begin
	 init_th(e1);
	 e0:=epot/tnodes;
	 e1:=e1/tnodes;
      end;
      if not(flip_debug) then begin
	 writeln('Final U=',e0:1:6,' U_test=',e1:1:6,' [eV/at]');
	 rewrite(txef,'Ef');
	 writeln(txef,e1*natt:1:5);
	 close(txef);
      end;
   end;
   close(txti);
   if deco then close(txta);
   close(recu);
   if save_xyzs then close(txyzs);
   if (deco) then begin
      rewrite(txt,'occupancy.out');
      for i:=1 to nnele do begin
	 writeln(txt,ele[i]:1);
	 for j:=1 to maxorb do begin
	    write(txt,j:2,'  ');
	    for k:=1 to loops do write(txt,o_ch[k,j,ele[i]]/o_al[k,j]:7:4);
	    writeln(txt);
	 end;
      end;
      writeln(txt,'0'); {vacancy}
      for j:=1 to maxorb do begin
	 write(txt,j:2,'  ');
	 for k:=1 to loops do write(txt,o_ch[k,j,0]/o_al[k,j]:7:4);
	 writeln(txt);
      end;
      close(txt);
   end;
end; {d04_main}
procedure density_composition;
var a : v33; det,vol:myreal;
   o  : integer;
begin
   if eval_nom and not(rchemx) then begin
      writeln('Density-composition: cannot evaluate without fractional chemistry.');
      halt;
   end;
   if rchemx and not(eval_nom) then begin
      writeln('Density-composition: cannot evaluate without atomic density.');
   end;
   if (rchemx and eval_nom) then begin
      dum:=0.0; for i:=1 to nnele do dum:=dum+xmele[i];
      if (abs(dum-1.0)>0.0001) then begin
	 writeln('Inconsistent composition in d04.inp: Sum(x)=',dum:1:4);
	 halt;
      end;
   end;
   A:=bl3;
   DET:=A[1,1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*A[3,1]
   +A[2,1]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*A[3,1]
   -A[1,1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*A[2,1];
   vol:=abs(det); {*multi}
   if debug[11] then writeln('UnitCell volume: ',vol:1:3);
   if not(cheminp) then begin
      if (tinp=0) then write('Decoration rule -> ') else
	 write('Config_para -> ');
      occupied_sites(natoc);
      for o:=1 to nnele do begin
	 mele_nom[o]:=mele_oc[o];
	 xmele[o]:=mele_oc[o]/natoc;
      end;
      natnom:=natoc;
      natlim:=natoc;
   end else begin
      write('chemistry redefinition -> ');
      if rchemx then begin
	 for o:=1 to nnele do begin
	    mele_nom[o]:=round(xmele[o]*rho*vol);
	 end;
	 natnom:=0; for o:=1 to nnele do natnom:=natnom+mele_nom[o];
      end else begin
	 for o:=1 to nnele do xmele[o]:=mele_nom[o]/natnom;
      end;
      rho_lim:=natnom/vol;
      natlim:=natnom;
   end;
   if (natnom=0) then begin
      write('Nominal number of atoms is 0. ');
      if cheminp then writeln('Chemistry defined by decoration rule.') else
	 writeln('Check chemistry in d04.inp.');
      halt;
   end;
   for i:=1 to nnele do pele[ele[i]]:=i;
   writeln('Rho_nom=',natnom/vol:1:5,', NAT_nom=',natnom:1,', ',nnele:1,'  elements:');
   for o:=1 to nnele do
      writeln(ele[o]:1,' ',xmele[o]:1:4,' ',mele_nom[o]:1);
end; { density_composition }
procedure startread(Var txt : text; Var x:char; Var eofi:boolean);
begin
   eofi:=false;
   repeat
      if (eoln(txt)) then readln(txt);
      eofi:=eof(txt);
   until not(eoln(txt)) or eofi;
   if not(eofi) then begin
      read(txt,x);
      if ((x<>'=') and (x<>'<') and (x<>'+')) then begin
	 readln(txt); eofi:=eof(txt);
	 if not(eofi) then startread(txt,x,eofi);
      end;
   end;
end; { startread }
procedure init_debug_options;
begin
   for i:=1 to 1000 do begin debug[i]:=false; debug_msg[i]:=' ??? '; end;
   {debug options}
   debug_msg[1]:='energy (flipping)';
   debug_msg[2]:='flipping insight';
   debug_msg[3]:='test update energy';
   debug_msg[4]:='atoms-init info + add_chem_para';
   debug_msg[5]:='swap (LG) energy update';
   debug_msg[6]:='g-flip energy update';
   debug_msg[7]:='tiling connectivity';
   debug_msg[8]:='relaxation';
   debug_msg[9]:='sym.-breaking';
   debug_msg[10]:='input files';
   debug_msg[11]:='INFO (general)';
   debug_msg[12]:='Export copa.ini';
   debug_msg[13]:='DFIT';
   debug_msg[22]:='MORE flipping info';
   debug_msg[23]:='explicit paths';
end;
procedure read_d04inp(d04inpfile :mystring);
var
   i,j,k,o : integer;
   txr	   : text;
begin
   hiter:=1;
   viter:=1;
   nwarn:=0;
   flip_debug:=false;
   flip_man:=false;
   flip_rec:=false;
   swap_man:=false;
   skip_reject:=false;
   table_reinit:=false;
   save_xyzs:=false;
   cheminp:=false;
   til1:=0; til2:=0;
   elimo:=0.0;
   erec1:=-100.0; erec2:=100.0;
   nodum:=false;
   durecalc:=false;
   export_err:=true;
   init_debug_options;
   level:=1; {surface}
   tlevel:=1;
   rcutn:=20.0;
   funi:=1; feps:=1e-6;
   nflipnu:=0;
   dfe:=true; {true: export delta_U type output for flip_replay}
   test_glob:=1; test_perm:=0; {wont use anymore}
   save_e:=10000; seed:=1;
   loops:=0;
   lgchp:=false; chvar:=false; natch:=0; chpot:=0.0;
   for o:=0 to 110 do begin chp[o]:=0; buflim[o]:=0; end; {chemical potentials; new}
   for o:=0 to nele do chpt[o]:=0.0; loops_ini:=0; {old}
   rchemi:=false; rchemx:=false; eval_nom:=false;
   for o:=1 to t_orb do o_frozen[o]:=false;
   for o:=1 to t_orb do o_rel[o]:=true; {relax status of an orbit}
   relax_man:=false;
   for i:=1 to tdim do
      for j:=1 to at_of_node do anb[i,j]:=false; {boolean for relax_man}
   nmrel:=0; nmrelf:=0;
   reset(txr,d04inpfile);
   readln(txr); {comment line}
   readln(txr,nschemes); {scheme}
   sname[1]:='LatticeGas';
   sname[2]:='SwapPair';
   sname[3]:='OrbitFlip';
   sname[4]:='ObjectFlip';
   sname[5]:='NodeFlip';
   sname[6]:='Swap2Pairs';
   sname[7]:='RelaxSymco';
   sname[8]:='DFIT';
   sname[9]:='CMC';
   sname[10]:='InstanceFlip';
   omcs:=0; tmcs:=0; smcs:=0; lmcs:=0; gmcs:=0; imcs:=0;
   nbond:=0; smcs2:=0; nbpa:=0; relmcs:=0;
   for i:=1 to nschemes do begin
      read(txr,sch[i,1],sch[i,2]);
      if (sch[i,1]=1) then lmcs:=lmcs+sch[i,2];
      if (sch[i,1]=2) then begin {swap}
	 nbond:=nbond+1;
	 smcs:=smcs+sch[i,2];
	 for o:=1 to 2 do read(txr,bondch[nbond,o]);
	 for o:=1 to 2 do read(txr,bondr[nbond,o]);
      end;
      if (sch[i,1]=3) then gmcs:=gmcs+sch[i,2];
      if (sch[i,1]=4) then omcs:=omcs+sch[i,2];
      if (sch[i,1]=5) then tmcs:=tmcs+sch[i,2];
      if (sch[i,1]=6) then begin
	 nbpa:=nbpa+1;
	 smcs2:=smcs2+sch[i,2];
	 for o:=1 to 2 do read(txr,twopa[nbpa,o]);
      end;
      if (sch[i,1]=7) then begin
	 read(txr,xgamma,dramax);
	 relmcs:=relmcs+sch[i,2];
      end;
      if (sch[i,1]=9) then begin
	 read(txr,sigmc);
	 relmcs:=relmcs+sch[i,2];
      end;
      if (sch[i,1]=10) then imcs:=imcs+sch[i,2];
      readln(txr);
   end;
   {check swap2 def is OK}
   for i:=1 to nbpa do 
      for o:=1 to 2 do if (twopa[i,o]>nbond) then begin
	 writeln('Inconsistent definition of swap2pairs: pair[',i:1,',',o:1,']>nbond, nbond=',nbond:1);
	 halt;
      end;
   readln(txr,viter,hiter); {MC_schedule}
   readln(txr,loops);
   for i:=1 to loops do readln(txr,mcs[i],beta[i]);
   readln(txr,o); {flip_record}
   if (o=1) then flip_rec:=true;
   read(txr,o); {flip_man}
   if (o>0) then begin
      flip_man:=true;
      nseq:=o;
      read(txr,o);
      nschemes:=1;
      sch[1,1]:=o;
      sch[1,2]:=1;
      for i:=1 to nseq do read(txr,vflip[i]);
      viter:=1; hiter:=1; loops:=1;
   end;
   readln(txr);
   read(txr,o); {swap_man}
   if (o>0) then begin
      swap_man:=true;
      nswapm:=o;
      for i:=1 to nswapm do read(txr,ato1[i],ato2[i]);
   end;
   readln(txr);
   read(txr,o); {relax_list}
   if (o>0) then begin
      read(txr,relmcs0,xgamma,dramax);
      relax_man:=true;
      nmrel:=o;
      for i:=1 to nmrel do read(txr,atr[i]); 
   end;
   readln(txr);
   read(txr,nflipnu); {flip_nu}
   for i:=1 to nflipnu do begin
      read(txr,flip_nu[i,1],flip_nu[i,2]);
   end;
   readln(txr);
   readln(txr,o); if(o=1) then sigmaflip:=true else sigmaflip:=false; {bowtie flip}
   {writeln('sigmaflip ',sigmaflip);}

   readln(txr,level,tlevel); {1=surface; n=add n boundaries around the surface}
   readln(txr,couple); {couple fat+skinny hex. flip}
   readln(txr,vpin); {pin a vertex}
   readln(txr,cueff); {ueff}
   readln(txr,rcutn); {rcutn}
   readln(txr,rcplus); {rcplus}
   readln(txr,feps); {flip_union_eps}
   readln(txr,seed); {random seed}
   randomize;
   readln(txr,o); {save_xyzs}
   if (o=1) then save_xyzs:=true;
   readln(txr,save_e); {save E every n-th step}
   readln(txr,hiterw);
   read(txr,o); {debug}
   for i:=1 to o do begin
      read(txr,j); debug[j]:=true;
   end;
   readln(txr);
   for i:=1 to 1000 do if debug[i] then writeln('Debug option ',i:1,': ',debug_msg[i]);
   readln(txr,o); {chemistry redef}
   if (o>0) then begin
      cheminp:=true;
      read(txr,m);
      if (m=0) then begin {integer-based chemistry}
	 rchemx:=false;
	 eval_nom:=false;
	 rchemi:=true;
	 read(txr,nnele);
	 for o:=1 to nnele do read(txr,ele[o],mele_nom[o]);
	 nat:=0; for o:=1 to nnele do nat:=nat+mele_nom[o];
	 for o:=1 to nnele do xmele[o]:=mele_nom[o]/nat;
	 {writeln('READINPDBG : ',mele_nom[1]:5,mele_nom[2]:5,nat:5);}
	 natnom:=nat;
	 natlim:=nat;
      end;
      if (m=1) then begin
	 {figure out integer representation after reading uc.inp.}
	 eval_nom:=true;
	 read(txr,rho,rho_lim);
      end;
      if (m=2) then begin
	 eval_nom:=true;
	 rchemx:=true;
	 read(txr,rho,nnele);
	 for i:=1 to nnele do read(txr,ele[i],xmele[i]);
      end;
      if (m=3) then begin
	 loops_ini:=loops_ini+1;
	 read(txr,mcs[loops+loops_ini],beta[loops+loops_ini]);
      end;
      if (m=4) then begin {old way of including chemical potentials}
	 lgchp:=true;
	 read(txr,o);
	 for i:=1 to o do read(txr,chpt[ele[i]]);
      end;
      readln(txr);
   end;
   read(txr,o); {chemical potentials}
   if (o>0) then begin
      chvar:=true;
      for m:=1 to o do begin
	 read(txr,i,cmp,k);
	 chp[i]:=cmp;
	 buflim[i]:=k; {buffer limit}
      end;
   end; readln(txr);
   readln(txr,o); if (o=1) then dodfit:=true else dodfit:=false;
   close(txr);
end; { read_d04inp }
procedure swap2(node1,atom1,node2,atom2	:integer );
Var i : integer;
begin
   if (node1>tnodes) or (node2>tnodes) then begin
      if (node1>tnodes) then writeln('SWAP_MAN: node1 > tnodes.');
      if (node2>tnodes) then writeln('SWAP_MAN: node2 > tnodes.');
      halt;
   end;
   if (atom1>nan[node1]) then begin
      writeln('SWAP_MAN: atom1 > nan[atom1].'); halt;
   end;
   if (atom2>nan[node2]) then begin
      writeln('SWAP_MAN: atom2 > nan[atom2].'); halt;
   end;
   if (n_a[node1,atom1,1]=n_a[node2,atom2,1]) then begin
      writeln('SWAP_MAN: chemistries of the two atoms are equal (',n_a[node1,atom1,1]:1,')'); halt;
   end;
   i:=n_a[node1,atom1,1];
   n_a[node1,atom1,1]:=n_a[node2,atom2,1];
   n_a[node2,atom2,1]:=i;
   para2xyz(n_para[node1,atom1],n_a[node1,atom1,1],1,1,n_r[node1,atom1]);
   para2xyz(n_para[node2,atom2],n_a[node2,atom2,1],1,1,n_r[node2,atom2]);
   writeln('SWAP_MAN: swapping ',node1:1,':',atom1:1,':',n_a[node1,atom1,1]:1,' <-> ',node2:1,':',atom2:1,':',n_a[node2,atom2,1]:1);
end; { swap2 }
BEGIN
   deconame:='';
   reset(txt,".d04inp");
   readln(txt,ticoname);
   readln(txt,tinp);
   readln(txt,o); read_tvl:=(o=1);
   readln(txt,inpname);
   close(txt);
   flip_replay:=(tinp=2);
   natch:=0;
   emin:=1.0e10;
   iname[1]:='ATOMIC DECORATION / PAIR POTENTIALS';
   iname[2]:='TILING OBJECTS / TILE HAMILTONIAN';
   BASIS; {define basis vectors}
   READ_D04INP(inpname);
   if debug[11] then begin
      writeln('DECOTYPE=',tinp:1);
   end;
   if (not(tinp=2) and flip_replay) then begin
      writeln('Flip-replay mode needs flip.rec input type.');
      writeln('This must be explicitly said in the 3rd line of the input file.');
      halt;
   end;
   i:=0; {addobj=0};
   READ_DECO(i); {reads objects.def, atoms.def and cond. orbit_labels.def}
   {READ_TILING}
   reset(txt,ticoname);
   IMPORT_TILING_BODY(tinp,txt);
   TABLE_INIT;
   lpar[1]:=1.0; lpar[2]:=multi; bl3[3,3]:=lpar[2];
   if deco then begin
      UNIT_CELL(multi); {reads uc.inp; MUST come after read_tiling}
      PP_INP; {reads pp.inp}
   end;
   for i:=1 to tnodes do nan[i]:=0;
   if (instances and not(read_tvl)) then begin
      for i:=1 to tnodes do new_pat_instance(i,0);
   end;
   flip_instances:=false;
   for i:=1 to 1000 do begin dbgs[i]:=debug[i]; debug[i]:=false; end;
   for i:=1 to tnodes do NEW_NODE_OBJECTS(i,1);
   for i:=1 to 1000 do debug[i]:=dbgs[i];
   if (tinp=1) then begin
      IMPORT_CONFIG(txt);
   end;
   if not(flip_replay) then close(txt);
   if deco then begin
      DENSITY_COMPOSITION;
   end;
   if (flip_replay) then begin
      if deco then begin
	 write('Flip-replay mode is designed to run in ');
	 writeln(iname[2],' mode.'); halt;
      end;
      readln(txt,couple);
      nseq:=0; repeat nseq:=nseq+1;
	 readln(txt,vflip[nseq],fene[nseq],dfene[nseq],doflip[nseq]);
      until eof(txt);
      close(txt);
      viter:=1; hiter:=1; loops:=1;
      nschemes:=1;
      sch[1,1]:=5; sch[1,2]:=1;
      fene_tot:=0.0;
      for o:=1 to n_obj do osta[o]:=0;
      rewrite(txfout,'flip.rep');
      writeln(txfout,n_obj:1,' 0  # of objects');
      writeln('FLIP-REPLAY: ',nseq:1,' flip(s) -> flip.rep.');
   end;
   table_reinit:=((gmcs>0) or (smcs>0) or (smcs2>0) or ((lmcs>0) and (lgchp)));
   export_occup:=(deco and not((tmcs>0) or swap_man or flip_man or flip_replay or relax_man));
   if export_occup then writeln('Exporting atomic occupations.');
   if debug[11] then begin
      writeln(ticoname,' : ',tnodes:1,' nodes');
      writeln('nu_offs=',nu_offs:1);
      writeln('Unit_cell_multi=',multi:1);
      if (multi>1) then writeln('Z-Tile-Hamiltonian enabled.');
      writeln('TMCS LMCS SMCS SMCS2 GMCS ',tmcs:1,' ',lmcs:1,' ',smcs:1,' ',smcs2:1,' ',gmcs:1);
      writeln('RCUT=',rcut:1:1,' [A]');
   end;
   if (tnodes=0) then begin
      writeln('0 tiling nodes. (Is MULTI defined in the input tiling?)');
      halt;
   end;
   export_obj_st('objects.ini'); {save init objects statistics}
   {NOTE: with swap_manual, d04_main is not RUN!}
   if dodfit then begin
      READ_DFIT_INP(dmcsteps,itmp0,io,jo); {only dmcsteps is taken the rest is ignored}
      if (dmcsteps>0)then begin
	 nschemes:=nschemes+1;
	 sch[nschemes,1]:=10;
	 sch[nschemes,2]:=dmcsteps;
      end;
      DFIT_INIT;
      nntable(rcplus);
      rewrite(txrepdf,'dfit.rep');
      rewrite(txpardf,'dfit.copas');
      export_tiling_header(txpardf);
      export_tiling_body(txpardf);
      writeln(txpardf,'start_paradump');
   end;
   if not(swap_man) then begin
      D04_MAIN;
      OCH_STATUS;
   end else begin {SWAP_MANUAL}
      occupied_sites(p);
      UPPN(e0);
      if(chvar) then begin
	 cmp:=0;
	 chempot(cmp);
	 e0:=e0+cmp;
      end;
      writeln('Initial U=',e0/p:1:6,' [eV/atom].');
      o:=0; for i:=1 to tnodes do
	 for j:=1 to nan[i] do begin o:=o+1;
	    for k:=1 to nswapm do begin
	       if (o=ato1[k]) then begin no1[k]:=i; at1[k]:=j; ato1[k]:=0; end;
	       if (o=ato2[k]) then begin no2[k]:=i; at2[k]:=j; ato2[k]:=0; end;
	    end;
	 end;
      for i:=1 to nswapm do begin
	 if (ato1[i]<>0) then begin
	    writeln('Could not find atom1 ',ato1[i]:1);
	    halt;
	 end;
	 if (ato2[i]<>0) then begin
	    writeln('Could not find atom2 ',ato2[i]:1);
	    halt;
	 end;
      end;
      for i:=1 to nswapm do swap2(no1[i],at1[i],no2[i],at2[i]);
      UPPN(e1);
      if(chvar) then begin
	 chempot(cmp);
	 e1:=e1+cmp;
      end;
      writeln('dU = ',(e1-e0):1:6,' [eV]. Final U=',e1/p:1:6,' [ev/atom].');
   end;
   if relax_man then begin {manual relaxation}
      occupied_sites(p);
      export_xyz('xyz.ini');
      UPPN(e0);
      o:=0;
      for i:=1 to tnodes do
	 for j:=1 to nan[i] do begin o:=o+1;
	    m:=0; for k:=1 to nmrel do if (o=atr[k]) then begin
	       anb[i,j]:=true; atr[k]:=0;
	    end;
	 end;
      for i:=1 to nmrel do if (atr[i]<>0) then begin
	 writeln('Could not find atom ',atr[i]:1,'.');
	 halt;
      end;
      if (nmrel=0) then begin
	 writeln('RELAX_MAN: no atoms on the input: ALL will be relaxed.');
	 for i:=1 to tnodes do
	    for j:=1 to nan[i] do begin nmrel:=nmrel+1; anb[i,j]:=true; end;
      end;
      writeln('RELAX_MAN: ',nmrel:1,' atoms to relax by ',relmcs0:1,' steps.');
      if (relmcs0>0) then begin
	 {relax(relmcs0,xgamma,dramax,epot);}
	 writeln('After ',relmcs0:1,' relaxation steps:');
	 writeln('dU = ',(epot-e0):1:6,' [eV]. Final relaxed U=',epot/p:1:6,' [ev/atom].');
	 {restore the unrelaxed state}
	 export_xyz('xyz.relm');
	 writeln('RELAX_MAN: exporting xyz.relm.');
	 for i:=1 to tnodes do
	    for j:=1 to nan[i] do if anb[i,j] then 
	       for o:=1 to 3 do n_r[i,j,o]:=n_r[i,j,o]-n_dr[i,j,o];
	 {export list of largest displacements}
	 writeln('RELAX_MAN: exporting drlist.out.');
	 rewrite(txt,'drlist.out');
	 writeln(txt,nmrel:1,'  relaxed atoms; # atom node node_atom chem orb |dR|[A] dZ[A]');
	 p:=0; for i:=1 to tnodes do
	    for j:=1 to nan[i] do begin
	       p:=p+1;
	       if anb[i,j] then begin
		  ortho3(n_dr[i,j],w3);
		  cmp:=absv3(w3);
		  writeln(txt,p:3,' ',i:3,' ',j:2,' ',n_a[i,j,1]:2,' ',n_a[i,j,2]:2,'   ',cmp:7:3,' ',w3[3]:7:3);
	       end;
	    end;
	 close(txt);
      end;
   end;
   if (dodfit) then begin
      close(txpardf);
      close(txrepdf);
   end;
   if deco then begin
      nntable(rcplus);
      cmp:=0.0;
      elastic_E(cmp,epot,elas,se,sel);
      export_xyz_enes('xyz.d04',se,sel);
   end;
   export_obj_st('objects.out'); {save final objects statistics}
   export_tiling('tiling.out');
   if (deco) then begin
      {if repre5 then export_atoms5d('atoms_5d.out');}
      {para_clean;}
      export_config('copa.out');
   end;
   if flip_rec then close(txfout);
END.

